Load the necessary libraries
library(car) #for regression diagnostics
library(broom) #for tidy output
library(ggfortify) #for model diagnostics
library(sjPlot) #for outputs
library(knitr) #for kable
library(effects) #for partial effects plots
library(ggeffects) #for partial effects plots
library(emmeans) #for estimating marginal means
library(MASS) #for glm.nb
library(MuMIn) #for AICc
library(tidyverse) #for data wrangling
library(nlme)
library(lme4)
library(glmmTMB)
library(broom.mixed)
library(glmmTMB) #for glmmTMB
library(DHARMa) #for residuals and diagnostics
library(performance) #for diagnostic plots
library(see) #for diagnostic plots
To investigate differential metabolic plasticity in barramundi (Lates calcarifer), Norin, Malte, and Clark (2015) exposed juvenile barramundi to various environmental changes (increased temperature, decreased salinity and increased hypoxia) as well as control conditions. Metabolic plasticity was calculated as the percentage difference in standard metabolic rate between the various treatment conditions and the standard metabolic rate under control conditions. They were interested in whether there was a relationship between metabolic plasticity and typical (control) metabolism and how the different treatment conditions impact on this relationship.
A total of 60 barramundi juveniles were subject to each of the three conditions (high temperature, low salinity and hypoxia) in addition to control conditions. Fish mass was also recorded as a covariate as this is known to influence metabolic parameters.
Barramundi
Sampling design
Format of norin.csv data files
| FISHID | MASS | TRIAL | SMR_contr | CHANGE |
|---|---|---|---|---|
| 1 | 35.69 | LowSalinity | 5.85 | -31.92 |
| 2 | 33.84 | LowSalinity | 6.53 | 2.52 |
| 3 | 37.78 | LowSalinity | 5.66 | -6.28 |
| .. | .. | .. | .. | .. |
| 1 | 36.80 | HighTemperature | 5.85 | 18.32 |
| 2 | 34.98 | HighTemperature | 6.53 | 19.06 |
| 3 | 38.38 | HighTemperature | 5.66 | 19.03 |
| .. | .. | .. | .. | .. |
| 1 | 45.06 | Hypoxia | 5.85 | -18.61 |
| 2 | 43.51 | Hypoxia | 6.53 | -5.37 |
| 3 | 45.11 | Hypoxia | 5.66 | -13.95 |
| FISHID | Categorical listing of the individual fish that are repeatedly sampled |
| MASS | Mass (g) of barramundi. Covariate in analysis |
| TRIAL | Categorical listing of the trial (LowSalinity: 10ppt salinity; HighTemperature: 35 degrees; Hypoxia: 45% air-sat. oxygen. |
| SMR_contr | Standard metabolic rate (mg/h/39.4 g of fish) under control trial conditions (35 ppt salinity, 29 degrees, normoxia) |
| CHANGE | Percentage difference in Standard metabolic rate (mg/h/39.4 g of fish) between Trial conditions and control adjusted for 'regression to the mean'. |
norin = read_csv('../public/data/norin.csv', trim_ws=TRUE)
glimpse(norin)
## Rows: 180
## Columns: 5
## $ FISHID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
## $ MASS <dbl> 35.69, 33.84, 37.78, 26.58, 37.62, 37.68, 30.62, 50.37, 24.…
## $ TRIAL <chr> "LowSalinity", "LowSalinity", "LowSalinity", "LowSalinity",…
## $ SMR_contr <dbl> 5.847466, 6.530707, 5.659556, 6.278200, 4.407336, 4.818589,…
## $ CHANGE <dbl> -31.919389, 2.520929, -6.284968, -4.346675, -3.071329, -15.…
Since the response variable is a difference between the SMR under control conditions compared to the specific trial conditions (either Low salinity, High temporature or Hypoxia), the only possible distribution we can use is a Gaussian. Well to clarify this, we could use a t-distribution, but this is not really available for frequentist mixed effects routines in R.
Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of temperature and (centered) mean fish size on SDA peak. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual fish.
Lets begin by declaring the categorical predictors and random effects as factors.
norin = norin %>% mutate(FISHID=factor(FISHID),
TRIAL=factor(TRIAL))
Each of the fish were exposed to each of the three trials (treatments). So lets explore the distributions of responses within each of these trials.
ggplot(norin, aes(y=CHANGE, x=TRIAL)) + geom_boxplot()
Conclusions:
The researchers considered that physiological plasticity might be effected by the fishes basal metabolic rate. For example, a fish with relatively high metabolism might have less scope to increase this metabolism than a fish with a lower metabolism. Therefore, the SMR under control conditions (just prior to the specific trial) could be added as a continuous covariate.
Doing so would introduce two additional assumptions:
ggplot(norin, aes(y=CHANGE, x=SMR_contr, shape=TRIAL, color=TRIAL)) +
geom_smooth(method='lm') + geom_point()
ggplot(norin, aes(y=CHANGE, x=SMR_contr, shape=TRIAL, color=TRIAL)) +
geom_smooth() + geom_point()
It might also be worth exploring how consistent the trial effect is across the different fish. This can give us an idea of whether the addition of a random intercept/slope model would be useful.
ggplot(norin, aes(y=CHANGE, x=as.numeric(FISHID), color=TRIAL)) +
geom_point() + geom_line()
Conclusions:
Finally, as an acknowledgement that the metabolic response might be influenced
by the mass of the fish, the researchers contemplated including fish Mass as a continuous covariate. There are numerous alternative ways to incorporate covariates that are known to impact on a response.
Lets explore the relationship between the response and Mass separately for each Trial.
ggplot(norin, aes(y=CHANGE, x=MASS, color=TRIAL)) +
geom_point() +
geom_smooth(method='lm')
Conclusions:
For each of the following modelling alternatives, we will:
##Compare models that estimate partial slope for MASS vs an offset for MASS
##must use ML to compare models that vary in fixed effects
norin.lme1 <- lme(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + offset(MASS), random=~1|FISHID, data=norin, method='ML')
## Unfortunately, update() does not remove offset().
## We will just have to write the other model out in full as well.
norin.lme2 <- lme(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + MASS, random=~1|FISHID, data=norin, method='ML')
## Now without MASS altogether
norin.lme3 <- lme(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE), random=~1|FISHID, data=norin, method='ML')
## Compare these models via AICc
AICc(norin.lme1,norin.lme2, norin.lme3)
## Alternatively, we can use sequential Likelihood Ratio Tests (LRT)
anova(norin.lme1, norin.lme2, norin.lme3)
Conclusions:
norin.lme3a = update(norin.lme3, method='REML')
norin.lme3b = update(norin.lme3a, random=~TRIAL|FISHID)
AICc(norin.lme3a, norin.lme3b)
anova(norin.lme3a, norin.lme3b)
Conclusions:
##Compare models that estimate partial slope for MASS vs an offset for MASS
##must use ML to compare models that vary in fixed effects
norin.lmer1 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + offset(MASS) + (1|FISHID), data=norin, REML=FALSE)
## Unfortunately, update() does not remove offset().
## We will just have to write the other model out in full as well.
norin.lmer2 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + MASS + (1|FISHID), data=norin, REML=FALSE)
## Now without MASS altogether
norin.lmer3 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE)+ (1|FISHID), data=norin, REML=FALSE)
## Compare these models via AICc
AICc(norin.lmer1,norin.lmer2, norin.lmer3)
## Alternatively, we can use sequential Likelihood Ratio Tests (LRT)
anova(norin.lmer1, norin.lmer2, norin.lmer3)
Conclusions:
norin.lmer3a = update(norin.lmer3, REML=TRUE)
## unfortunately, the following random intercept/slope model cannot be fit
##norin.lmer3b = update(norin.lmer3a, ~ . -(1|FISHID) +(TRIAL|FISHID))
##AICc(norin.lmer3a, norin.lmer3b)
##anova(norin.lmer3a, norin.lmer3b)
Conclusions:
##Compare models that estimate partial slope for MASS vs an offset for MASS
norin.glmmTMB1 = glmmTMB(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + offset(MASS) + (1|FISHID), data=norin, REML=FALSE)
## Unfortunately, update() does not remove offset().
## We will just have to write the other model out in full as well.
norin.glmmTMB2 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + MASS + (1|FISHID), data=norin, REML=FALSE)
## Now without MASS altogether
norin.glmmTMB3 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + (1|FISHID), data=norin, REML=FALSE)
## Compare these models via AICc
AICc(norin.glmmTMB1,norin.glmmTMB2, norin.glmmTMB3)
Conclusions:
norin.glmmTMB3a <- update(norin.glmmTMB3, REML=TRUE)
## unfortunately, the following random intercept/slope model cannot be fit
## norin.glmmTMB3b <- update(norin.glmmTMB3a, ~ . - (1|FISHID)+ (TRIAL|FISHID))
##AICc(norin.glmmTMB1a, norin.glmmTMB1b)
Conclusions:
Conclusions:
plot_grid(plot_model(norin.lme3b, type='diag')[-2])
Conclusions:
performance::check_model(norin.lme3b)
## Could not compute standard errors from random effects for diagnostic plot.
## Homogeneity of variance could not be computed. Cannot extract residual variance from objects of class 'lme'.
#norin.resid = simulateResiduals(norin.lme3a, plot=TRUE)
plot_grid(plot_model(norin.lmer3a, type='diag')[-2])
Conclusions:
performance::check_model(norin.lmer3a)
Conclusions:
norin.resid = simulateResiduals(norin.lmer3a, plot=TRUE)
Conclusions:
plot_grid(plot_model(norin.glmmTMB3a, type='diag')[-2])
Conclusions:
performance::check_model(norin.glmmTMB3a)
Conclusions:
norin.resid = simulateResiduals(norin.glmmTMB3a, plot=TRUE)
Conclusions:
plot_model(norin.lme3b, type='eff', terms=c('SMR_contr', 'TRIAL'), show.data=TRUE)
We can also explore the fixed effects (as a caterpillar plot).
plot_model(norin.lme3a, type='est')
And similarly, the individual random effects (intercepts of each fish).
#plot_model(norin.lme3a, type='re')
plot(allEffects(norin.lme3b), multiline=TRUE, ci.style='bands')
ggpredict(norin.lme3b, c('SMR_contr', 'TRIAL')) %>% plot(add.data=TRUE)
#ggemmeans(norin.lme3a, ~SMR_contr*TRIAL) %>% plot(add.data=TRUE)
plot_model(norin.lmer3a, type='eff', show.data=TRUE) %>% plot_grid
plot_model(norin.lmer3a, type='eff', terms=c('SMR_contr', 'TRIAL'), show.data=TRUE)
We can also explore the fixed effects (as a caterpillar plot).
plot_model(norin.lmer3a, type='est')
And similarly, the individual random effects (intercepts of each fish).
plot_model(norin.lmer3a, type='re')
plot(allEffects(norin.lmer3a), multiline=TRUE, ci.style='bands')
ggpredict(norin.lmer3a, c('SMR_contr', 'TRIAL')) %>% plot(add.data=TRUE)
ggemmeans(norin.lmer3a, ~SMR_contr*TRIAL) %>% plot(add.data=TRUE)
plot_model(norin.glmmTMB3a, type='eff', show.data=TRUE) %>% plot_grid
plot_model(norin.glmmTMB3a, type='eff', terms=c('SMR_contr', 'TRIAL'), show.data=TRUE)
We can also explore the fixed effects (as a caterpillar plot).
plot_model(norin.glmmTMB3a, type='est')
And similarly, the individual random effects (intercepts of each fish).
plot_model(norin.glmmTMB3a, type='re')
plot(allEffects(norin.glmmTMB3a), multiline=TRUE, ci.style='bands')
ggpredict(norin.glmmTMB3a, c('SMR_contr', 'TRIAL')) %>% plot(add.data=TRUE)
ggemmeans(norin.glmmTMB3a, ~SMR_contr*TRIAL) %>% plot(add.data=TRUE)
summary(norin.lme3b)
## Linear mixed-effects model fit by REML
## Data: norin
## AIC BIC logLik
## 1593.961 1635.028 -783.9803
##
## Random effects:
## Formula: ~TRIAL | FISHID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 16.323848 (Intr) TRIALH
## TRIALHypoxia 14.926864 -0.601
## TRIALLowSalinity 30.304028 -0.111 -0.215
## Residual 8.659992
##
## Fixed effects: CHANGE ~ TRIAL * scale(SMR_contr, scale = FALSE)
## Value Std.Error DF
## (Intercept) 50.38244 2.385594 116
## TRIALHypoxia -50.22621 2.492663 116
## TRIALLowSalinity -42.22318 4.219647 116
## scale(SMR_contr, scale = FALSE) -21.55293 3.820913 58
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 13.51107 3.992402 116
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 10.18925 6.758444 116
## t-value p-value
## (Intercept) 21.119453 0.0000
## TRIALHypoxia -20.149618 0.0000
## TRIALLowSalinity -10.006330 0.0000
## scale(SMR_contr, scale = FALSE) -5.640779 0.0000
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 3.384195 0.0010
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 1.507633 0.1344
## Correlation:
## (Intr) TRIALH TRIALL s(Ss=F
## TRIALHypoxia -0.620
## TRIALLowSalinity -0.215 -0.035
## scale(SMR_contr, scale = FALSE) 0.000 0.000 0.000
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 0.000 0.000 0.000 -0.620
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 0.000 0.000 0.000 -0.215
## TRIALHs=F
## TRIALHypoxia
## TRIALLowSalinity
## scale(SMR_contr, scale = FALSE)
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -0.035
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.66585731 -0.30692682 0.02265975 0.29410064 1.72690964
##
## Number of Observations: 180
## Number of Groups: 60
Conclusions:
intervals(norin.lme3a)
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est.
## (Intercept) 44.1837220 50.38244
## TRIALHypoxia -58.2125517 -50.22621
## TRIALLowSalinity -50.2095159 -42.22318
## scale(SMR_contr, scale = FALSE) -31.5868953 -21.55293
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 0.7196578 13.51107
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -2.6021570 10.18925
## upper
## (Intercept) 56.58115
## TRIALHypoxia -42.23987
## TRIALLowSalinity -34.23684
## scale(SMR_contr, scale = FALSE) -11.51896
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 26.30248
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 22.98066
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: FISHID
## lower est. upper
## sd((Intercept)) 5.931668 9.996282 16.84613
##
## Within-group standard error:
## lower est. upper
## 19.41879 22.08543 25.11827
tidy(norin.lme3a, conf.int=TRUE)
tidy(norin.lme3a, conf.int=TRUE) %>% kable
| effect | group | term | estimate | std.error | df | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | fixed | (Intercept) | 50.382438 | 3.129675 | 116 | 16.098295 | 0.0000000 | 44.1837220 | 56.58115 |
| fixed | fixed | TRIALHypoxia | -50.226212 | 4.032230 | 116 | -12.456187 | 0.0000000 | -58.2125517 | -42.23987 |
| fixed | fixed | TRIALLowSalinity | -42.223176 | 4.032230 | 116 | -10.471420 | 0.0000000 | -50.2095159 | -34.23684 |
| fixed | fixed | scale(SMR_contr, scale = FALSE) | -21.552927 | 5.012679 | 58 | -4.299682 | 0.0000663 | -31.5868953 | -11.51896 |
| fixed | fixed | TRIALHypoxia:scale(SMR_contr, scale = FALSE) | 13.511068 | 6.458266 | 116 | 2.092058 | 0.0386147 | 0.7196578 | 26.30248 |
| fixed | fixed | TRIALLowSalinity:scale(SMR_contr, scale = FALSE) | 10.189253 | 6.458266 | 116 | 1.577707 | 0.1173561 | -2.6021570 | 22.98066 |
| ran_pars | FISHID | sd_(Intercept) | 9.996282 | NA | NA | NA | NA | 5.9316684 | 16.84613 |
| ran_pars | Residual | sd_Observation | 22.085434 | NA | NA | NA | NA | NA | NA |
Conclusions:
# warning this is only appropriate for html output
sjPlot::tab_model(norin.lme3a, show.se=TRUE, show.aic=TRUE)
| CHANGE | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 50.38 | 3.13 | 44.18 – 56.58 | <0.001 |
| TRIAL [Hypoxia] | -50.23 | 4.03 | -58.21 – -42.24 | <0.001 |
| TRIAL [LowSalinity] | -42.22 | 4.03 | -50.21 – -34.24 | <0.001 |
| SMR_contr | -21.55 | 5.01 | -31.59 – -11.52 | <0.001 |
|
TRIAL [Hypoxia] * SMR_contr |
13.51 | 6.46 | 0.72 – 26.30 | 0.039 |
|
TRIAL [LowSalinity] * SMR_contr |
10.19 | 6.46 | -2.60 – 22.98 | 0.117 |
| Random Effects | ||||
| σ2 | 487.77 | |||
| τ00 FISHID | 99.93 | |||
| N FISHID | 60 | |||
| Observations | 180 | |||
| Marginal R2 / Conditional R2 | 0.541 / NA | |||
| AIC | 1636.349 | |||
summary(norin.lmer3a)
## Linear mixed model fit by REML ['lmerMod']
## Formula: CHANGE ~ TRIAL * scale(SMR_contr, scale = FALSE) + (1 | FISHID)
## Data: norin
##
## REML criterion at convergence: 1620.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8551 -0.6207 0.0224 0.5421 3.7421
##
## Random effects:
## Groups Name Variance Std.Dev.
## FISHID (Intercept) 99.93 9.996
## Residual 487.77 22.085
## Number of obs: 180, groups: FISHID, 60
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 50.382 3.130 16.098
## TRIALHypoxia -50.226 4.032 -12.456
## TRIALLowSalinity -42.223 4.032 -10.471
## scale(SMR_contr, scale = FALSE) -21.553 5.013 -4.300
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 13.511 6.458 2.092
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 10.189 6.458 1.578
##
## Correlation of Fixed Effects:
## (Intr) TRIALH TRIALL s(Ss=F TRIALHs=F
## TRIALHypoxi -0.644
## TRIALLwSlnt -0.644 0.500
## s(SMR_,s=FA 0.000 0.000 0.000
## TRIALH:(s=F 0.000 0.000 0.000 -0.644
## TRIALLS:s=F 0.000 0.000 0.000 -0.644 0.500
Conclusions:
confint(norin.lmer3a)
## 2.5 % 97.5 %
## .sig01 3.1895683 14.79565
## .sigma 19.2319058 24.78039
## (Intercept) 44.3172310 56.44765
## TRIALHypoxia -58.0590014 -42.39342
## TRIALLowSalinity -50.0559656 -34.39039
## scale(SMR_contr, scale = FALSE) -31.2673343 -11.83852
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 0.9655934 26.05654
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -2.3562215 22.73473
tidy(norin.lmer3a, conf.int=TRUE)
tidy(norin.lmer3a, conf.int=TRUE) %>% kable
| effect | group | term | estimate | std.error | statistic | conf.low | conf.high |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 50.382438 | 3.129675 | 16.098295 | 44.2483875 | 56.51649 |
| fixed | NA | TRIALHypoxia | -50.226212 | 4.032230 | -12.456186 | -58.1292378 | -42.32319 |
| fixed | NA | TRIALLowSalinity | -42.223176 | 4.032230 | -10.471420 | -50.1262020 | -34.32015 |
| fixed | NA | scale(SMR_contr, scale = FALSE) | -21.552927 | 5.012679 | -4.299682 | -31.3775983 | -11.72826 |
| fixed | NA | TRIALHypoxia:scale(SMR_contr, scale = FALSE) | 13.511068 | 6.458267 | 2.092058 | 0.8530984 | 26.16904 |
| fixed | NA | TRIALLowSalinity:scale(SMR_contr, scale = FALSE) | 10.189253 | 6.458267 | 1.577707 | -2.4687165 | 22.84722 |
| ran_pars | FISHID | sd__(Intercept) | 9.996281 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 22.085435 | NA | NA | NA | NA |
Conclusions:
# warning this is only appropriate for html output
sjPlot::tab_model(norin.lmer3a, show.se=TRUE, show.aic=TRUE)
| CHANGE | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 50.38 | 3.13 | 44.25 – 56.52 | <0.001 |
| TRIAL [Hypoxia] | -50.23 | 4.03 | -58.13 – -42.32 | <0.001 |
| TRIAL [LowSalinity] | -42.22 | 4.03 | -50.13 – -34.32 | <0.001 |
| SMR_contr | -21.55 | 5.01 | -31.38 – -11.73 | <0.001 |
|
TRIAL [Hypoxia] * SMR_contr |
13.51 | 6.46 | 0.85 – 26.17 | 0.036 |
|
TRIAL [LowSalinity] * SMR_contr |
10.19 | 6.46 | -2.47 – 22.85 | 0.115 |
| Random Effects | ||||
| σ2 | 487.77 | |||
| τ00 FISHID | 99.93 | |||
| ICC | 0.17 | |||
| N FISHID | 60 | |||
| Observations | 180 | |||
| Marginal R2 / Conditional R2 | 0.494 / 0.580 | |||
| AIC | 1636.349 | |||
summary(norin.glmmTMB3a)
## Linear mixed model fit by REML ['lmerMod']
## Formula: CHANGE ~ TRIAL * scale(SMR_contr, scale = FALSE) + (1 | FISHID)
## Data: norin
##
## REML criterion at convergence: 1620.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8551 -0.6207 0.0224 0.5421 3.7421
##
## Random effects:
## Groups Name Variance Std.Dev.
## FISHID (Intercept) 99.93 9.996
## Residual 487.77 22.085
## Number of obs: 180, groups: FISHID, 60
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 50.382 3.130 16.098
## TRIALHypoxia -50.226 4.032 -12.456
## TRIALLowSalinity -42.223 4.032 -10.471
## scale(SMR_contr, scale = FALSE) -21.553 5.013 -4.300
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 13.511 6.458 2.092
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 10.189 6.458 1.578
##
## Correlation of Fixed Effects:
## (Intr) TRIALH TRIALL s(Ss=F TRIALHs=F
## TRIALHypoxi -0.644
## TRIALLwSlnt -0.644 0.500
## s(SMR_,s=FA 0.000 0.000 0.000
## TRIALH:(s=F 0.000 0.000 0.000 -0.644
## TRIALLS:s=F 0.000 0.000 0.000 -0.644 0.500
Conclusions:
confint(norin.glmmTMB3a)
## 2.5 % 97.5 %
## .sig01 3.1895683 14.79565
## .sigma 19.2319058 24.78039
## (Intercept) 44.3172310 56.44765
## TRIALHypoxia -58.0590014 -42.39342
## TRIALLowSalinity -50.0559656 -34.39039
## scale(SMR_contr, scale = FALSE) -31.2673343 -11.83852
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 0.9655934 26.05654
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -2.3562215 22.73473
tidy(norin.glmmTMB3a, conf.int=TRUE)
tidy(norin.glmmTMB3a, conf.int=TRUE) %>% kable
| effect | group | term | estimate | std.error | statistic | conf.low | conf.high |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 50.382438 | 3.129675 | 16.098295 | 44.2483875 | 56.51649 |
| fixed | NA | TRIALHypoxia | -50.226212 | 4.032230 | -12.456186 | -58.1292378 | -42.32319 |
| fixed | NA | TRIALLowSalinity | -42.223176 | 4.032230 | -10.471420 | -50.1262020 | -34.32015 |
| fixed | NA | scale(SMR_contr, scale = FALSE) | -21.552927 | 5.012679 | -4.299682 | -31.3775983 | -11.72826 |
| fixed | NA | TRIALHypoxia:scale(SMR_contr, scale = FALSE) | 13.511068 | 6.458267 | 2.092058 | 0.8530984 | 26.16904 |
| fixed | NA | TRIALLowSalinity:scale(SMR_contr, scale = FALSE) | 10.189253 | 6.458267 | 1.577707 | -2.4687165 | 22.84722 |
| ran_pars | FISHID | sd__(Intercept) | 9.996281 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 22.085435 | NA | NA | NA | NA |
Conclusions:
# warning this is only appropriate for html output
sjPlot::tab_model(norin.glmmTMB3a, show.se=TRUE, show.aic=TRUE)
| CHANGE | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 50.38 | 3.13 | 44.25 – 56.52 | <0.001 |
| TRIAL [Hypoxia] | -50.23 | 4.03 | -58.13 – -42.32 | <0.001 |
| TRIAL [LowSalinity] | -42.22 | 4.03 | -50.13 – -34.32 | <0.001 |
| SMR_contr | -21.55 | 5.01 | -31.38 – -11.73 | <0.001 |
|
TRIAL [Hypoxia] * SMR_contr |
13.51 | 6.46 | 0.85 – 26.17 | 0.036 |
|
TRIAL [LowSalinity] * SMR_contr |
10.19 | 6.46 | -2.47 – 22.85 | 0.115 |
| Random Effects | ||||
| σ2 | 487.77 | |||
| τ00 FISHID | 99.93 | |||
| ICC | 0.17 | |||
| N FISHID | 60 | |||
| Observations | 180 | |||
| Marginal R2 / Conditional R2 | 0.494 / 0.580 | |||
| AIC | 1636.349 | |||
The presence of a significant interaction suggests that the nature of the relationship between SMR change and the control SMR differs between the three trials (High Temperature, Hypoxia and Low Salinity). Equivalently, the diffences between the three trials is not constant over all levels of control SMR.
To further tease these interactions apart, we could either:
We will now explore each of the above, although typically we would only do one or the other (depending on which variable was of more interest).
emtrends(norin.lme3a, pairwise~TRIAL, var='SMR_contr')
## $emtrends
## TRIAL SMR_contr.trend SE df lower.CL upper.CL
## HighTemperature -21.55 5.01 58 -31.6 -11.52
## Hypoxia -8.04 5.01 58 -18.1 1.99
## LowSalinity -11.36 5.01 58 -21.4 -1.33
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia -13.51 6.46 116 -2.092 0.0959
## HighTemperature - LowSalinity -10.19 6.46 116 -1.578 0.2594
## Hypoxia - LowSalinity 3.32 6.46 116 0.514 0.8645
##
## Degrees-of-freedom method: containment
## P value adjustment: tukey method for comparing a family of 3 estimates
Conclusions:
We will estimated the pairwise differences between each of the three trials at the minimum, mean and maximum control SMR.
norin.grid <- with(norin, list(SMR_contr=c(min(SMR_contr), mean(SMR_contr), max(SMR_contr))))
norin.grid
## $SMR_contr
## [1] 3.984319 5.178857 6.700537
emmeans(norin.lme3a, pairwise~TRIAL|SMR_contr, at=norin.grid)
## $emmeans
## SMR_contr = 3.98:
## TRIAL emmean SE df lower.CL upper.CL
## HighTemperature 76.128 6.76 58 62.60 89.65
## Hypoxia 9.763 6.76 58 -3.76 23.29
## LowSalinity 21.734 6.76 58 8.21 35.26
##
## SMR_contr = 5.18:
## TRIAL emmean SE df lower.CL upper.CL
## HighTemperature 50.382 3.13 59 44.12 56.64
## Hypoxia 0.156 3.13 59 -6.11 6.42
## LowSalinity 8.159 3.13 59 1.90 14.42
##
## SMR_contr = 6.70:
## TRIAL emmean SE df lower.CL upper.CL
## HighTemperature 17.586 8.24 58 1.08 34.09
## Hypoxia -12.081 8.24 58 -28.58 4.42
## LowSalinity -9.133 8.24 58 -25.64 7.37
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## SMR_contr = 3.98:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 66.37 8.70 116 7.624 <.0001
## HighTemperature - LowSalinity 54.39 8.70 116 6.249 <.0001
## Hypoxia - LowSalinity -11.97 8.70 116 -1.375 0.3572
##
## SMR_contr = 5.18:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 50.23 4.03 116 12.456 <.0001
## HighTemperature - LowSalinity 42.22 4.03 116 10.471 <.0001
## Hypoxia - LowSalinity -8.00 4.03 116 -1.985 0.1205
##
## SMR_contr = 6.70:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 29.67 10.62 116 2.793 0.0167
## HighTemperature - LowSalinity 26.72 10.62 116 2.515 0.0351
## Hypoxia - LowSalinity -2.95 10.62 116 -0.278 0.9584
##
## Degrees-of-freedom method: containment
## P value adjustment: tukey method for comparing a family of 3 estimates
Conclusions:
r.squaredGLMM(norin.lme3a)
## R2m R2c
## [1,] 0.4942091 0.5802091
## Nakagawa's R2
performance::r2_nakagawa(norin.lme3a)
## Random effect variances not available. Returned R2 does not account for random effects.
## # R2 for Mixed Models
##
## Conditional R2: NA
## Marginal R2: 0.541
emtrends(norin.lmer3a, pairwise~TRIAL, var='SMR_contr')
## $emtrends
## TRIAL SMR_contr.trend SE df lower.CL upper.CL
## HighTemperature -21.55 5.01 164 -31.5 -11.66
## Hypoxia -8.04 5.01 164 -17.9 1.86
## LowSalinity -11.36 5.01 164 -21.3 -1.47
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia -13.51 6.46 116 -2.092 0.0959
## HighTemperature - LowSalinity -10.19 6.46 116 -1.578 0.2594
## Hypoxia - LowSalinity 3.32 6.46 116 0.514 0.8645
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
norin.emt <- emtrends(norin.lmer3a, pairwise~TRIAL, var='SMR_contr')$contrasts %>% as.data.frame
Conclusions:
norin.grid <- with(norin, list(SMR_contr=c(min(SMR_contr), mean(SMR_contr), max(SMR_contr))))
norin.grid
## $SMR_contr
## [1] 3.984319 5.178857 6.700537
emmeans(norin.lmer3a, pairwise~TRIAL|SMR_contr, at=norin.grid)
## $emmeans
## SMR_contr = 3.98:
## TRIAL emmean SE df lower.CL upper.CL
## HighTemperature 76.128 6.76 164 62.79 89.47
## Hypoxia 9.763 6.76 164 -3.58 23.10
## LowSalinity 21.734 6.76 164 8.39 35.07
##
## SMR_contr = 5.18:
## TRIAL emmean SE df lower.CL upper.CL
## HighTemperature 50.382 3.13 164 44.20 56.56
## Hypoxia 0.156 3.13 164 -6.02 6.34
## LowSalinity 8.159 3.13 164 1.98 14.34
##
## SMR_contr = 6.70:
## TRIAL emmean SE df lower.CL upper.CL
## HighTemperature 17.586 8.24 164 1.31 33.87
## Hypoxia -12.081 8.24 164 -28.36 4.20
## LowSalinity -9.133 8.24 164 -25.41 7.15
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## SMR_contr = 3.98:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 66.37 8.70 116 7.624 <.0001
## HighTemperature - LowSalinity 54.39 8.70 116 6.249 <.0001
## Hypoxia - LowSalinity -11.97 8.70 116 -1.375 0.3572
##
## SMR_contr = 5.18:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 50.23 4.03 116 12.456 <.0001
## HighTemperature - LowSalinity 42.22 4.03 116 10.471 <.0001
## Hypoxia - LowSalinity -8.00 4.03 116 -1.985 0.1205
##
## SMR_contr = 6.70:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 29.67 10.62 116 2.793 0.0167
## HighTemperature - LowSalinity 26.72 10.62 116 2.515 0.0351
## Hypoxia - LowSalinity -2.95 10.62 116 -0.278 0.9584
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
Conclusions:
r.squaredGLMM(norin.lmer3a)
## R2m R2c
## [1,] 0.4942091 0.580209
## Nakagawa's R2
performance::r2_nakagawa(norin.lmer3a)
## # R2 for Mixed Models
##
## Conditional R2: 0.580
## Marginal R2: 0.494
emtrends(norin.glmmTMB3a, pairwise~TRIAL, var='SMR_contr')
## $emtrends
## TRIAL SMR_contr.trend SE df lower.CL upper.CL
## HighTemperature -21.55 5.01 164 -31.5 -11.66
## Hypoxia -8.04 5.01 164 -17.9 1.86
## LowSalinity -11.36 5.01 164 -21.3 -1.47
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia -13.51 6.46 116 -2.092 0.0959
## HighTemperature - LowSalinity -10.19 6.46 116 -1.578 0.2594
## Hypoxia - LowSalinity 3.32 6.46 116 0.514 0.8645
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
norin.emt <- emtrends(norin.glmmTMB3a, pairwise~TRIAL, var='SMR_contr')$contrasts %>% as.data.frame
Conclusions:
norin.grid <- with(norin, list(SMR_contr=c(min(SMR_contr), mean(SMR_contr), max(SMR_contr))))
norin.grid
## $SMR_contr
## [1] 3.984319 5.178857 6.700537
emmeans(norin.glmmTMB3a, pairwise~TRIAL|SMR_contr, at=norin.grid)
## $emmeans
## SMR_contr = 3.98:
## TRIAL emmean SE df lower.CL upper.CL
## HighTemperature 76.128 6.76 164 62.79 89.47
## Hypoxia 9.763 6.76 164 -3.58 23.10
## LowSalinity 21.734 6.76 164 8.39 35.07
##
## SMR_contr = 5.18:
## TRIAL emmean SE df lower.CL upper.CL
## HighTemperature 50.382 3.13 164 44.20 56.56
## Hypoxia 0.156 3.13 164 -6.02 6.34
## LowSalinity 8.159 3.13 164 1.98 14.34
##
## SMR_contr = 6.70:
## TRIAL emmean SE df lower.CL upper.CL
## HighTemperature 17.586 8.24 164 1.31 33.87
## Hypoxia -12.081 8.24 164 -28.36 4.20
## LowSalinity -9.133 8.24 164 -25.41 7.15
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## SMR_contr = 3.98:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 66.37 8.70 116 7.624 <.0001
## HighTemperature - LowSalinity 54.39 8.70 116 6.249 <.0001
## Hypoxia - LowSalinity -11.97 8.70 116 -1.375 0.3572
##
## SMR_contr = 5.18:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 50.23 4.03 116 12.456 <.0001
## HighTemperature - LowSalinity 42.22 4.03 116 10.471 <.0001
## Hypoxia - LowSalinity -8.00 4.03 116 -1.985 0.1205
##
## SMR_contr = 6.70:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 29.67 10.62 116 2.793 0.0167
## HighTemperature - LowSalinity 26.72 10.62 116 2.515 0.0351
## Hypoxia - LowSalinity -2.95 10.62 116 -0.278 0.9584
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
Conclusions:
r.squaredGLMM(norin.glmmTMB3a)
## R2m R2c
## [1,] 0.4942091 0.580209
## Nakagawa's R2
performance::r2_nakagawa(norin.glmmTMB3a)
## # R2 for Mixed Models
##
## Conditional R2: 0.580
## Marginal R2: 0.494
norin.grid=with(norin, list(SMR_contr=seq(min(SMR_contr), max(SMR_contr), len=100)))
newdata = emmeans(norin.lme3a, ~SMR_contr|TRIAL, at=norin.grid) %>%
as.data.frame
head(newdata)
ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
geom_ribbon(aes(ymin=lower.CL, ymax=upper.CL, fill=TRIAL), alpha=0.3) +
geom_line(aes(, color=TRIAL)) +
theme_classic() +
theme(legend.position = c(0.99, 0.99),
legend.justification = c(1, 1))
obs <- norin.lme3a %>%
augment() %>%
mutate(PartialObs=.fitted + .resid)
ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
geom_ribbon(aes(ymin=lower.CL, ymax=upper.CL, fill=TRIAL), alpha=0.3) +
geom_line(aes(, color=TRIAL)) +
geom_point(data=obs, aes(y=PartialObs, color=TRIAL)) +
theme_classic()
norin.grid=with(norin, list(SMR_contr=seq(min(SMR_contr), max(SMR_contr), len=100)))
newdata = emmeans(norin.lmer3a, ~SMR_contr|TRIAL, at=norin.grid) %>%
as.data.frame
head(newdata)
ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
geom_ribbon(aes(ymin=lower.CL, ymax=upper.CL, fill=TRIAL), alpha=0.3) +
geom_line(aes(, color=TRIAL)) +
theme_classic() +
theme(legend.position = c(0.99, 0.99),
legend.justification = c(1, 1))
obs <- norin.lmer3a %>%
augment() %>%
bind_cols(norin %>% select(SMR_contr)) %>%
mutate(PartialObs=.fitted + .resid)
ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
geom_ribbon(aes(ymin=lower.CL, ymax=upper.CL, fill=TRIAL), alpha=0.3) +
geom_line(aes(, color=TRIAL)) +
geom_point(data=obs, aes(y=PartialObs, color=TRIAL)) +
theme_classic()
norin.grid=with(norin, list(SMR_contr=seq(min(SMR_contr), max(SMR_contr), len=100)))
newdata = emmeans(norin.glmmTMB3a, ~SMR_contr|TRIAL, at=norin.grid) %>%
as.data.frame
head(newdata)
ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
geom_ribbon(aes(ymin=lower.CL, ymax=upper.CL, fill=TRIAL), alpha=0.3) +
geom_line(aes(, color=TRIAL)) +
theme_classic() +
theme(legend.position = c(0.99, 0.99),
legend.justification = c(1, 1))
obs <- norin.glmmTMB3a %>%
augment() %>%
bind_cols(norin %>% select(SMR_contr)) %>%
mutate(PartialObs=.fitted + .resid)
ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
geom_ribbon(aes(ymin=lower.CL, ymax=upper.CL, fill=TRIAL), alpha=0.3) +
geom_line(aes(, color=TRIAL)) +
geom_point(data=obs, aes(y=PartialObs, color=TRIAL)) +
theme_classic()
Norin, Tommy, Hans Malte, and Timothy D. Clark. 2015. “Differential Plasticity of Metabolic Rate Phenotypes in a Tropical Fish Facing Environmental Change.” Functional Ecology 30 (3): 369–78. https://doi.org/10.1111/1365-2435.12503.